#include "bootstrapping.h"
Bootstrapper::Bootstrapper(Context &context, Scheme &scheme, SecretKey &secretKey1, EvalModPoly &evalModPoly, EncodingMatrix &encodingMatrix)
        : context(context), scheme(scheme), secretKey(secretKey1), evalModPoly(evalModPoly), encodingMatrix(encodingMatrix)
{
    Q = 0x4001b00001;
    levelStart = 16;
    scalingFactor = 1L << 36;
    sineType = SineType::Cos2;
    messageRatio = 4;
    k = 10;
    sineDeg = 31;
    doubleAngle = 2;
    arcSineDeg = 0;
    // --------初始化evalmod中间参数---------
    Ciphertext ct0;
    for (int i = 0; i < int(log2(sineDeg + 1)); i++) {
        evalSinePolyPool.push_back(new Ciphertext());
        evalBasis.push_back(ct0);
    }
    for (int i = 0; i < sineDeg + 1; i++) {
        chebyshevPolyCoeffTreePool.push_back(new Polynomial());
        chebyshevTreePool.push_back(new Ciphertext());
    }
    // --------初始化evalmod中间参数---------
}

/*
 * @brief SecretKey 类的构造函数，用于初始化一个新的密钥对象。
 *
 * @param context：加密方案的上下文对象，包含了加密过程中需要的各种参数。
 *
 */
void Bootstrapper::SplitCoeffsPolyVector(int split, vector<Polynomial *> chebyshev_coeff_tree_pool, int tree_idx)
{
    Polynomial *coeffs = chebyshev_coeff_tree_pool[tree_idx];
    Polynomial *coeffsq = chebyshev_coeff_tree_pool[tree_idx * 2];
    Polynomial *coeffsr = chebyshev_coeff_tree_pool[tree_idx * 2 + 1];

    coeffsr->Coeffs = vector<complex<double>>(split);
    if (coeffs->MaxDeg == coeffs->Degree()) {
        coeffsr->MaxDeg = split - 1;
    } else {
        coeffsr->MaxDeg = coeffs->MaxDeg - (coeffs->Degree() - split + 1);
    }

    for (int i = 0; i < split; i++)
        coeffsr->Coeffs[i] = coeffs->Coeffs[i];

    coeffsq->Coeffs = vector<complex<double>>(coeffs->Degree() - split + 1);
    coeffsq->MaxDeg = coeffs->MaxDeg;

    coeffsq->Coeffs[0] = coeffs->Coeffs[split];

    for (int i = split + 1, j = 1; i < coeffs->Degree() + 1; i++, j++) {
        coeffsq->Coeffs[i - split] = complex<double>(coeffs->Coeffs[i].real() * 2, 0);
        coeffsr->Coeffs[split - j] = coeffsr->Coeffs[split - j] - coeffs->Coeffs[i];
    }
}

/*
 * @brief ：递归准备Chebyshev系数树，根据给定的拆分日志和多项式度数进行计算
 *
 * @param logSplit ：拆分的日志值
 * @param logDegree ：多项式的日志度数
 * @param tree_idx ：树的索引，用于在系数树池中定位相应的多项式
 */
void Bootstrapper::PrepareChebyshevCoeffsTree(int logSplit, int logDegree, int tree_idx)
{
    Polynomial *poly = chebyshevPolyCoeffTreePool[tree_idx];
    if (poly->Degree() < (1 << logSplit)) {
        if (logSplit > 1 && poly->MaxDeg % (1 << (logSplit + 1)) > (1 << (logSplit - 1))) {
            logDegree = ceil(log2(poly->Degree()));
            logSplit = logDegree >> 1;

            PrepareChebyshevCoeffsTree(logSplit, logDegree, tree_idx);
            return;
        }
        return;
    }
    int nextPower = 1 << logSplit;
    for (nextPower; nextPower < (poly->Degree() >> 1) + 1;)
        nextPower <<= 1;

    Polynomial *coeffsq = chebyshevPolyCoeffTreePool[tree_idx * 2];
    Polynomial *coeffsr = chebyshevPolyCoeffTreePool[tree_idx * 2 + 1];

    SplitCoeffsPolyVector(nextPower, chebyshevPolyCoeffTreePool, tree_idx);

    PrepareChebyshevCoeffsTree(logSplit, logDegree, tree_idx * 2);
    PrepareChebyshevCoeffsTree(logSplit, logDegree, tree_idx * 2 + 1);
}

/*
 * @brief ：准备正弦Chebyshev多项式，使用预先计算的系数生成系数树
 */
void Bootstrapper::PrepareSineChebyshevPoly()
{
    Polynomial eval_sine_chebyshev(evalSineChebyshevCoeff);
    chebyshevPolyCoeffTreePool[1]->Coeffs = eval_sine_chebyshev.Coeffs;
    chebyshevPolyCoeffTreePool[1]->MaxDeg = eval_sine_chebyshev.MaxDeg;

    int degree = chebyshevPolyCoeffTreePool[1]->Degree();
    int logDegree = ceil(log2(degree));
    int logSplit = (logDegree >> 1);
    PrepareChebyshevCoeffsTree(logSplit, logDegree, 1);
}

/*
 * @brief ：准备正弦Chebyshev多项式评估，使用手动计算的系数生成系数树
 *
 * @param evalModPolyPre ：用于计算正弦Chebyshev多项式的评估多项式
 */
void Bootstrapper::PrepareSineChebyshevPolyEval(EvalModPoly evalModPolyPre)
{
    evalModPolyPre.sinePoly->Coeffs.emplace_back(0, 0);
    evalModPolyPre.sinePoly->MaxDeg = evalModPolyPre.sinePoly->MaxDeg + 1;
    chebyshevPolyCoeffTreePool[1]->Coeffs = evalModPolyPre.sinePoly->Coeffs;
    chebyshevPolyCoeffTreePool[1]->MaxDeg = evalModPolyPre.sinePoly->MaxDeg;

    int degree = chebyshevPolyCoeffTreePool[1]->Degree();
    int logDegree = ceil(log2(degree));
    int logSplit = (logDegree >> 1);
    PrepareChebyshevCoeffsTree(logSplit, logDegree, 1);
}

/*
 * @brief ：计算Chebyshev基，通过解密输入密文并进行多次操作，计算Chebyshev基的值
 *
 * @param ct0 ：输入密文
 * @param target_scale ：目标缩放因子
 * @param secretKeyCom ：私钥，用于解密
 */
void Bootstrapper::ComputeChebyshevBasis(Ciphertext &ct0, SecretKey &secretKeyCom)
{
    *(evalSinePolyPool[0]) = ct0;

    complex<double> *dvecbasic;
    complex<double> *dvecbasic1;
    complex<double> *mvecbasic1_pingfang;
    complex<double> *mvecbasic1_pingfang_2;
    complex<double> *mvecbasic1_pingfang_2_1;
    Ciphertext guofu;
    dvecbasic1 = scheme.Decrypt(secretKeyCom, ct0);
    for (int idx = 1; idx < ceil(log2(sineDeg)); idx++) {
        dvecbasic = scheme.Decrypt(secretKeyCom, *evalSinePolyPool[idx - 1]);

        dvecbasic1 = scheme.Decrypt(secretKeyCom, *evalSinePolyPool[idx - 1]);

        *(evalSinePolyPool[idx]) = scheme.Square(*(evalSinePolyPool[idx - 1]));
        scheme.ReScaleByAndEqual(*(evalSinePolyPool[idx]), 1);

        *(evalSinePolyPool[idx]) = scheme.Add(*(evalSinePolyPool[idx]), *(evalSinePolyPool[idx]));


        scheme.AddConstAndEqual(*(evalSinePolyPool[idx]), -1);
    }
}

/*
 * @brief ：在Chebyshev系数树上执行迭代计算，基于目标缩放因子和给定的私钥进行一系列计算
 *
 * @param target_scale ：目标缩放因子
 * @param logDegreeEval ：多项式度数的日志值
 * @param secretKeyEval ：私钥，用于解密
 */
void Bootstrapper::EvalIteration(__uint128_t target_scale, SecretKey &secretKeyEval)
{
    uint64_t currentQi = context.qVec[evalSinePolyPool[0]->l - 1];
    __uint128_t ctScale = target_scale * currentQi;
    __uint128_t constScale = ctScale / evalSinePolyPool[0]->scale;

    complex<double> *dvecbasic;
    complex<double> *dvecbasic1;
    for (int tree_idx = 31; tree_idx >= 16; tree_idx--) {
        Polynomial *poly = chebyshevPolyCoeffTreePool[tree_idx];

        double c1 = poly->Coeffs[0].real();
        double c2 = poly->Coeffs[1].real();


        chebyshevTreePool[tree_idx]->l = evalSinePolyPool[0]->l;
        chebyshevTreePool[tree_idx]->scale = evalSinePolyPool[0]->scale;

        dvecbasic = scheme.Decrypt(secretKeyEval, *(evalSinePolyPool[0])); // 解密结果

        *chebyshevTreePool[tree_idx] = scheme.MultByConst(*(evalSinePolyPool[0]), c2); // 正常的标量乘

        scheme.ReScaleByAndEqual(*chebyshevTreePool[tree_idx], 1);

        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx]); // 解密结果
        scheme.AddConstAndEqual(*chebyshevTreePool[tree_idx], c1);
        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx]); // 解密结果
        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx]); // 解密结果
    }

    for (int tree_idx = 15; tree_idx >= 1; tree_idx--) {
        int idx_nextPower = 4 - int(log2(tree_idx));
        Ciphertext *xpow = evalSinePolyPool[idx_nextPower];

        Polynomial *poly = chebyshevPolyCoeffTreePool[tree_idx];
        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx * 2]);      // 解密结果
        dvecbasic1 = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx * 2 + 1]); // 解密结果
        dvecbasic1 = scheme.Decrypt(secretKeyEval, *xpow);                                // 解密结果

        scheme.MultAndEqual(*chebyshevTreePool[tree_idx * 2], *xpow);
        scheme.ReScaleByAndEqual(*chebyshevTreePool[tree_idx * 2], 1);

        scheme.ModDownByAndEqual(*chebyshevTreePool[tree_idx * 2 + 1], 1);
        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx * 2]);      // 解密结果
        dvecbasic1 = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx * 2 + 1]); // 解密结果
        *chebyshevTreePool[tree_idx] =
                scheme.Add(*chebyshevTreePool[tree_idx * 2 + 1], *chebyshevTreePool[tree_idx * 2]);
        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx]); // 解密结果
        dvecbasic = scheme.Decrypt(secretKeyEval, *chebyshevTreePool[tree_idx]); // 解密结果
    }
}

/*
 * @brief ：执行同态取模操作，调整密文的层级并执行必要的后处理
 *
 * @param ct0 ：输入密文
 */
void Bootstrapper::EvalMod(Ciphertext &ct0)
{
    if (ct0.l > levelStart) {
        scheme.ModDownToAndEqual(ct0, levelStart);
    } else if (ct0.l < levelStart) {
        cout << " ！！密文层级小于evalmod开始层级！！ " << endl;
    }

    __uint128_t prevScaleCt = ct0.scale;

    //    ct0.scale = scalingFactor;

    __uint128_t targetScale = ct0.scale;
    //    for (int i = 0; i < evalModPoly.doubleAngle; i++) {
    //        targetScale = sqrt(targetScale *context.qVec[(evalModPoly.levelStart - evalModPoly.sinePoly->Depth() -
    //        evalModPoly.doubleAngle + i)]);
    //    }

    // Division by 1/2^r and change of variable for the Chebyshev evaluation  x-0.025
    scheme.AddConstAndEqual(ct0, -0.5 / (evalModPoly.scFac * (evalModPoly.sinePoly->B - evalModPoly.sinePoly->A)));
    // ------------------------------------------------------------

    /* ***********************chebyshev coeff tree************************** */
    PrepareSineChebyshevPoly(); // 使用提前存好的系数
    //    bootstrapper.PrepareSineChebyshevPolyEval(evalModPoly);//使用计算的系数

    /* ************************计算chebyshev基******************************* */
    ComputeChebyshevBasis(ct0, secretKey);

    /* ************************迭代计算树******************************* */
    EvalIteration(targetScale, secretKey);

    /* ************************结果存回ct0中******************************* */
    ct0 = *chebyshevTreePool[1];

    // -------------------evalmod后处理----------------------------
    double sqrt2Pi_temp = double(evalModPoly.sqrt2Pi);
    cout << "sqrt2Pi_temp: " << sqrt2Pi_temp << endl;

    for (int i = 0; i < evalModPoly.doubleAngle; i++) {
        sqrt2Pi_temp *= sqrt2Pi_temp;
        scheme.SquareAndEqual(ct0);
        scheme.ReScaleByAndEqual(ct0, 1);
        scheme.AddAndEqual(ct0, ct0);
        scheme.AddConstAndEqual(ct0, -sqrt2Pi_temp);
    }
}

enum { UPPER_BOUND_PS = 2204 };
const std::complex<double> I(0.0, 1.0);
double PREC = std::pow(2, -20);
inline bool IsNotEqualOne(double val)
{
    if (1 - PREC >= val) {
        return true;
    }
    if (1 + PREC <= val) {
        return true;
    }
    return false;
}

/*
 * @brief ：根据给定的输入密文和系数，递归地应用Paterson-Stockmeyer算法进行Chebyshev多项式的评估
 *
 * @param x ：输入密文
 * @param coefficients ：Chebyshev多项式的系数
 * @param kInner ：多项式的阶数
 * @param m ：多项式的层数
 * @param T ：存储临时结果的密文向量
 * @param T2 ：存储最终结果的密文向量
 * @return Ciphertext ：返回经过评估的密文
 */
Ciphertext Bootstrapper::InnerEvalChebyshevPS(Ciphertext x, const vector<double> &coefficients, uint32_t kInner,
                                              uint32_t m, vector<Ciphertext> &T, vector<Ciphertext> &T2)
{
    // Compute kInner*2^{m-1}-kInner because we use it a lot
    uint32_t k2m2k = kInner * (1 << (m - 1)) - kInner;

    // Divide coefficients by T^{kInner*2^{m-1}}
    std::vector<double> Tkm(int32_t(k2m2k + kInner) + 1, 0.0);
    Tkm.back() = 1;
    auto divqr = LongDivisionChebyshev(coefficients, Tkm);

    // Subtract x^{kInner(2^{m-1} - 1)} from r
    std::vector<double> r2 = divqr.r;
    if (int32_t(k2m2k - Degree(divqr.r)) <= 0) {
        r2[int32_t(k2m2k)] -= 1;
        r2.resize(Degree(r2) + 1);
    } else {
        r2.resize(int32_t(k2m2k + 1), 0.0);
        r2.back() = -1;
    }

    // Divide r2 by q
    auto divcs = LongDivisionChebyshev(r2, divqr.q);

    // Add x^{kInner(2^{m-1} - 1)} to s
    std::vector<double> s2 = divcs.r;
    s2.resize(int32_t(k2m2k + 1), 0.0);
    s2.back() = 1;

    // Evaluate c at u
    Ciphertext cu;
    uint32_t dc = Degree(divcs.q);
    bool flag_c = false;
    if (dc >= 1) {
        if (dc == 1) {
            if (divcs.q[1] != 1) {
                cu = scheme.MultByConst(T.front(), divcs.q[1]);
                scheme.ReScaleByAndEqual(cu, 1);
            } else {
                cu = T.front();
            }
        } else {
            std::vector<Ciphertext> ctxs(dc);
            std::vector<double> weights(dc);

            for (uint32_t i = 0; i < dc; i++) {
                ctxs[i] = T[i];
                weights[i] = divcs.q[i + 1];
            }

            cu = EvalLinearWSumMutable(ctxs, weights);
        }

        // adds the free term (at x^0)
        scheme.AddConstAndEqual(cu, divcs.q.front() / 2);
        // Need to reduce levels up to the level of T2[m-1].
        long levelDiff = T2[m - 1].l - cu.l;
        scheme.ModDownByAndEqual(cu, levelDiff);
        flag_c = true;
    }

    // Evaluate q and s2 at u. If their degrees are larger than kInner, then recursively apply the Paterson-Stockmeyer
    // algorithm.
    Ciphertext qu;

    if (Degree(divqr.q) > kInner) {
        qu = InnerEvalChebyshevPS(x, divqr.q, kInner, m - 1, T, T2);
    } else {
        // dq = kInner from construction
        // perform scalar multiplication for all other terms and sum them up if there are non-zero coefficients
        auto qcopy = divqr.q;
        qcopy.resize(kInner);
        if (Degree(qcopy) > 0) {
            std::vector<Ciphertext> ctxs(Degree(qcopy));
            std::vector<double> weights(Degree(qcopy));

            for (uint32_t i = 0; i < Degree(qcopy); i++) {
                ctxs[i] = T[i];
                weights[i] = divqr.q[i + 1];
            }

            qu = EvalLinearWSumMutable(ctxs, weights);
            // the highest order coefficient will always be a power of two up to 2^{m-1} because q is "monic" but the
            // Chebyshev rule adds a factor of 2 we don't need to increase the depth by multiplying the highest order
            // coefficient, but instead checking and summing, since we work with m <= 4.
            Ciphertext sum = T[kInner - 1];
            for (uint32_t i = 0; i < log2(divqr.q.back()); i++) {
                sum = scheme.Add(sum, sum);
            }

            AdjustLevelsAndDepthInPlace(qu, sum);
            scheme.AddAndEqual(qu, sum);
        } else {
            Ciphertext sum = T[kInner - 1];
            for (uint32_t i = 0; i < log2(divqr.q.back()); i++) {
                sum = scheme.Add(sum, sum);
            }
            qu = sum;
        }

        // adds the free term (at x^0)

        scheme.AddConstAndEqual(qu, divqr.q.front() / 2);
        // The number of levels of qu is the same as the number of levels of T[kInner-1] or T[kInner-1] + 1.
        // No need to reduce it to T2[m-1] because it only reaches here when m = 2.
    }

    Ciphertext su;

    if (Degree(s2) > kInner) {
        su = InnerEvalChebyshevPS(x, s2, kInner, m - 1, T, T2);
    } else {
        // ds = kInner from construction
        // perform scalar multiplication for all other terms and sum them up if there are non-zero coefficients
        auto scopy = s2;
        scopy.resize(kInner);
        if (Degree(scopy) > 0) {
            std::vector<Ciphertext> ctxs(Degree(scopy));
            std::vector<double> weights(Degree(scopy));

            for (uint32_t i = 0; i < Degree(scopy); i++) {
                ctxs[i] = T[i];
                weights[i] = s2[i + 1];
            }

            su = EvalLinearWSumMutable(ctxs, weights);
            // the highest order coefficient will always be 1 because s2 is monic.

            Ciphertext test;
            test = T[kInner - 1];
            AdjustLevelsAndDepthInPlace(su, test);
            scheme.AddAndEqual(su, test);
        } else {
            su = T[kInner - 1];
        }

        // adds the free term (at x^0)

        scheme.AddConstAndEqual(su, s2.front() / 2);
        // The number of levels of su is the same as the number of levels of T[kInner-1] or T[kInner-1] + 1. Need to
        // reduce it to
        scheme.ModDownByAndEqual(su, 1);
    }

    Ciphertext result;

    if (flag_c) {
        result = scheme.Add(T2[m - 1], cu);
    } else {
        result = scheme.AddConst(T2[m - 1], divcs.q.front() / 2);
    }


    scheme.MultAndEqual(result, qu);
    scheme.ReScaleByAndEqual(result, 1);

    AdjustLevelsAndDepthInPlace(result, su);

    scheme.AddAndEqual(result, su);

    return result;
}

/*
 * @brief 根据给定的输入密文和系数，递归地应用Paterson-Stockmeyer算法进行Chebyshev多项式的评估
 *
 * @param x 输入密文，表示多项式的自变量
 * @param coefficients Chebyshev多项式的系数
 * @param a 和 b 用于线性变换的常数
 * @return Ciphertext 返回经过评估后的密文
 */
Ciphertext Bootstrapper::EvalChebyshevSeriesPS(Ciphertext x, const vector<double> &coefficients, double a, double b)
{
    
    uint32_t n = Degree(coefficients);

    std::vector<double> f2 = coefficients;

    x.scale = context.p;

    // Make sure the coefficients do not have the zero dominant terms
    if (coefficients[coefficients.size() - 1] == 0)
        f2.resize(n + 1);

    std::vector<uint32_t> degs = ComputeDegreesPS(n);
    uint32_t k1 = degs[0];
    uint32_t m = degs[1];


    // computes linear transformation y = -1 + 2 (x-a)/(b-a)
    // consumes one level when a <> -1 && b <> 1

    std::vector<Ciphertext> T(k1);
    if ((a - std::round(a) < 1e-10) && (b - std::round(b) < 1e-10) && (std::round(a) == -1) && (std::round(b) == 1)) {
        // no linear transformation is needed if a = -1, b = 1

        T[0] = x;
    } else {
        // linear transformation is needed
        double alpha = 2 / (b - a);
        double beta = 2 * a / (b - a);


        T[0] = scheme.MultByConst(x, alpha);
        scheme.ReScaleByAndEqual(T[0], 1);
        scheme.AddConstAndEqual(T[0], -1.0 - beta);
    }

    Ciphertext y = T[0];

    // Computes Chebyshev polynomials up to degree k1
    // for y: T_1(y) = y, T_2(y), ... , T_k(y)
    // uses binary tree multiplication
    for (uint32_t i = 2; i <= k1; i++) {
        // if i is a power of two
        if (!(i & (i - 1))) {
            // compute T_{2i}(y) = 2*T_i(y)^2 - 1

            auto square = scheme.Square(T[i / 2 - 1]);
            scheme.ReScaleByAndEqual(square, 1);
            T[i - 1] = scheme.Add(square, square);
            scheme.AddConstAndEqual(T[i - 1], -1.0);
        } else {
            // non-power of 2
            if (i % 2 == 1) {
                // if i is odd
                // compute T_{2i+1}(y) = 2*T_i(y)*T_{i+1}(y) - y

                Ciphertext test1, test2;
                Ciphertext prod;
                prod = T[i / 2 - 1];
                test1 = T[i / 2 - 1];
                test2 = T[i / 2];
                AdjustLevelsAndDepthInPlace(prod, test2); // 归一输入层数
                scheme.MultAndEqual(prod, test2);

                T[i - 1] = scheme.Add(prod, prod);
                scheme.ReScaleByAndEqual(T[i - 1], 1);
                test1 = T[i - 1];
                test2 = y;
                AdjustLevelsAndDepthInPlace(T[i - 1], y); // 归一输入层数
                T[i - 1] = scheme.Sub(T[i - 1], y);
            } else {
                // i is even but not power of 2
                // compute T_{2i}(y) = 2*T_i(y)^2 - 1


                auto square = scheme.Square(T[i / 2 - 1]);
                T[i - 1] = scheme.Add(square, square);
                scheme.ReScaleByAndEqual(T[i - 1], 1);
                scheme.AddConstAndEqual(T[i - 1], -1.0);
            }
        }
    }
    int wssb = 1;

    if (wssb == 0) {
        // brings all powers of x to the same level
        for (size_t i = 1; i < k1; i++) {
            long levelDiff = T[k1 - 1].l - T[i - 1].l;
            scheme.ModDownByAndEqual(T[i - 1], levelDiff);
        }
    } else {
        for (size_t i = 1; i < k1; i++) {
            AdjustLevelsAndDepthInPlace(T[i - 1], T[k1 - 1]);
        }
    }

    std::vector<Ciphertext> T2(m);
    // Compute the Chebyshev polynomials T_k(y), T_{2k}(y), T_{4k}(y), ... , T_{2^{m-1}k1}(y)
    // T2[0] is used as a placeholder
    T2.front() = T.back();
    for (uint32_t i = 1; i < m; i++) {
        auto square = scheme.Square(T2[i - 1]);
        scheme.ReScaleByAndEqual(square, 1);
        T2[i] = scheme.Add(square, square);
        scheme.AddConstAndEqual(T2[i], -1.0);
    }

    // computes T_{k1(2*m - 1)}(y)
    auto T2km1 = T2.front();
    for (uint32_t i = 1; i < m; i++) {
        // compute T_{k1(2*m - 1)} = 2*T_{k1(2^{m-1}-1)}(y)*T_{k1*2^{m-1}}(y) - T_k(y)

        Ciphertext midTest;
        midTest = T2[i];
        AdjustLevelsAndDepthInPlace(T2km1, midTest); // 归一输入层数
        auto prod = T2km1;
        scheme.MultAndEqual(prod, midTest);
        T2km1 = scheme.Add(prod, prod);
        scheme.ReScaleByAndEqual(T2km1, 1);
        midTest = T2.front();
        AdjustLevelsAndDepthInPlace(T2km1, midTest); // 归一输入层数
        T2km1 = scheme.Sub(T2km1, midTest);
    }

    // We also need to reduce the number of levels of T[k1-1] and of T2[0] by another level.
    //  cc->LevelReduceInPlace(T[k1-1], nullptr);
    //  cc->LevelReduceInPlace(T2.front(), nullptr);

    // Compute k1*2^{m-1}-k1 because we use it a lot
    uint32_t k2m2k = k1 * (1 << (m - 1)) - k1;

    // Add T^{k1(2^m - 1)}(y) to the polynomial that has to be evaluated
    f2.resize(2 * k2m2k + k1 + 1, 0.0);
    f2.back() = 1;

    // Divide f2 by T^{k1*2^{m-1}}
    std::vector<double> Tkm(int32_t(k2m2k + k1) + 1, 0.0);
    Tkm.back() = 1;
    auto divqr = LongDivisionChebyshev(f2, Tkm);

    // Subtract x^{k1(2^{m-1} - 1)} from r
    std::vector<double> r2 = divqr.r;
    if (int32_t(k2m2k - Degree(divqr.r)) <= 0) {
        r2[int32_t(k2m2k)] -= 1;
        r2.resize(Degree(r2) + 1);
    } else {
        r2.resize(int32_t(k2m2k + 1), 0.0);
        r2.back() = -1;
    }

    // Divide r2 by q
    auto divcs = LongDivisionChebyshev(r2, divqr.q);

    // Add x^{k1(2^{m-1} - 1)} to s
    std::vector<double> s2 = divcs.r;
    s2.resize(int32_t(k2m2k + 1), 0.0);
    s2.back() = 1;

    // Evaluate c at u
    Ciphertext cu;
    uint32_t dc = Degree(divcs.q);
    bool flag_c = false;
    if (dc >= 1) {
        if (dc == 1) {
            if (divcs.q[1] != 1) {
                cu = scheme.MultByConst(T.front(), divcs.q[1]);
                scheme.ReScaleByAndEqual(cu, 1);
            } else {
                cu = T.front();
            }
        } else {
            std::vector<Ciphertext> ctxs(dc);
            std::vector<double> weights(dc);

            for (uint32_t i = 0; i < dc; i++) {
                ctxs[i] = T[i];
                weights[i] = divcs.q[i + 1];
            }

            cu = EvalLinearWSumMutable(ctxs, weights);
        }

        // adds the free term (at x^0)

        cu = scheme.AddConst(cu, divcs.q.front() / 2);
        // TODO : Andrey why not T2[m-1]->GetLevel() instead?
        // Need to reduce levels to the level of T2[m-1].


        flag_c = true;
    }

    // Evaluate q and s2 at u. If their degrees are larger than k1, then recursively apply the Paterson-Stockmeyer
    // algorithm.
    Ciphertext qu;

    if (Degree(divqr.q) > k1) {
        qu = InnerEvalChebyshevPS(x, divqr.q, k1, m - 1, T, T2);
    } else {
        // dq = k1 from construction
        // perform scalar multiplication for all other terms and sum them up if there are non-zero coefficients
        auto qcopy = divqr.q;
        qcopy.resize(k1);
        if (Degree(qcopy) > 0) {
            std::vector<Ciphertext> ctxs(Degree(qcopy));
            std::vector<double> weights(Degree(qcopy));

            for (uint32_t i = 0; i < Degree(qcopy); i++) {
                ctxs[i] = T[i];
                weights[i] = divqr.q[i + 1];
            }

            qu = EvalLinearWSumMutable(ctxs, weights);
            // the highest order coefficient will always be 2 after one division because of the Chebyshev division rule

            Ciphertext sum = scheme.Add(T[k1 - 1], T[k1 - 1]);
            AdjustLevelsAndDepthInPlace(qu, sum);
            scheme.AddAndEqual(qu, sum);
        } else {
            qu = T[k1 - 1];

            for (uint32_t i = 1; i < divqr.q.back(); i++) {
                Ciphertext test;
                test = T[k1 - 1];
                AdjustLevelsAndDepthInPlace(qu, test);

                scheme.AddAndEqual(qu, test);
            }
        }

        // adds the free term (at x^0)

        scheme.AddConstAndEqual(qu, divqr.q.front() / 2);
        // The number of levels of qu is the same as the number of levels of T[k1-1] + 1.
        // Will only get here when m = 2, so the number of levels of qu and T2[m-1] will be the same.
    }

    Ciphertext su;

    if (Degree(s2) > k1) {
        su = InnerEvalChebyshevPS(x, s2, k1, m - 1, T, T2);
    } else {
        // ds = k1 from construction
        // perform scalar multiplication for all other terms and sum them up if there are non-zero coefficients
        auto scopy = s2;
        scopy.resize(k1);
        if (Degree(scopy) > 0) {
            std::vector<Ciphertext> ctxs(Degree(scopy));
            std::vector<double> weights(Degree(scopy));

            for (uint32_t i = 0; i < Degree(scopy); i++) {
                ctxs[i] = T[i];
                weights[i] = s2[i + 1];
            }

            su = EvalLinearWSumMutable(ctxs, weights);
            // the highest order coefficient will always be 1 because s2 is monic.

            AdjustLevelsAndDepthInPlace(su, T[k1 - 1]);

            scheme.AddAndEqual(su, T[k1 - 1]);
        } else {
            su = T[k1 - 1];
        }

        // adds the free term (at x^0)

        scheme.AddConstAndEqual(su, s2.front() / 2);
        // The number of levels of su is the same as the number of levels of T[k1-1] + 1.
        // Will only get here when m = 2, so need to reduce the number of levels by 1.
    }


    Ciphertext result;

    if (flag_c) {
        result = scheme.Add(T2[m - 1], cu);
    } else {
        result = scheme.AddConst(T2[m - 1], divcs.q.front() / 2);
    }

    scheme.MultAndEqual(result, qu);

    scheme.ReScaleByAndEqual(result, 1);

    AdjustLevelsAndDepthInPlace(su, result);

    scheme.AddAndEqual(result, su);

    result = scheme.Sub(result, T2km1);


    return result;
}

/*
 * @brief 计算给定多项式系数的阶数
 *
 * @param coefficients 输入的多项式系数
 * @return uint32_t 返回多项式的阶数
 */
uint32_t Bootstrapper::Degree(const vector<double> &coefficients)
{
    const size_t coefficientsSize = coefficients.size();
    if (!coefficientsSize) {
        cout << ("The coefficients vector can not be empty") << endl;
    }

    int32_t indx = coefficientsSize;
    while (--indx >= 0) {
        if (coefficients[indx])
            break;
    }

    // indx becomes negative (-1) only when all coefficients are zeroes. in this case we return 0
    return static_cast<uint32_t>((indx < 0) ? 0 : indx);
}

/*
 * @brief 计算用于Paterson-Stockmeyer算法的多项式的k和m值
 *
 * @param n 多项式的阶数
 * @return std::vector<uint32_t> 返回一个包含k和m值的向量
 */
std::vector<uint32_t> Bootstrapper::ComputeDegreesPS(const uint32_t n)
{
    if (n == 0) {
        cout << "ComputeDegreesPS: The degree is zero. There is no need to evaluate the polynomial." << endl;
    }
    // index n-1 in the vector corresponds to degree n
    if (n <= UPPER_BOUND_PS) { // hard-coded values
        static const std::vector<uint32_t> mlist = PopulateParameterPS(UPPER_BOUND_PS);
        uint32_t m = mlist[n - 1];
        uint32_t k1 = std::floor(n / ((1 << m) - 1)) + 1;

        return std::vector<uint32_t>{ k1, m };
    } else { // heuristic for larger degrees
        std::vector<uint32_t> klist;
        std::vector<uint32_t> mlist;
        std::vector<uint32_t> multlist;

        for (uint32_t i = 1; i <= n; i++) {
            for (uint32_t m = 1; m <= std::ceil(log2(n / i) + 1) + 1; m++) {
                if (int32_t(n) - int32_t(i * ((1 << m) - 1)) < 0) {
                    if (std::abs(std::floor(log2(i)) - std::floor(log2(sqrt(n / 2)))) <= 1) {
                        klist.push_back(i);
                        mlist.push_back(m);
                        multlist.push_back(i + 2 * m + (1 << (m - 1)) - 4);
                    }
                }
            }
        }
        uint32_t minIndex = std::min_element(multlist.begin(), multlist.end()) - multlist.begin();

        return std::vector<uint32_t>{ klist[minIndex], mlist[minIndex] };
    }
}

/*
 * @brief 填充用于Paterson-Stockmeyer算法的参数m值
 *
 * @param upperBoundDegree 最高阶数
 * @return std::vector<uint32_t> 返回包含每个阶数对应的m值的向量
 */
std::vector<uint32_t> Bootstrapper::PopulateParameterPS(const uint32_t upperBoundDegree)
{
    std::vector<uint32_t> mlist(upperBoundDegree);

    std::fill(mlist.begin(), mlist.begin() + 2, 1);           // n in [1,2], m = 1
    std::fill(mlist.begin() + 2, mlist.begin() + 11, 2);      // n in [3,11], m = 2
    std::fill(mlist.begin() + 11, mlist.begin() + 13, 3);     // n in [12,13], m = 3
    std::fill(mlist.begin() + 13, mlist.begin() + 17, 2);     // n in [14,17], m = 2
    std::fill(mlist.begin() + 17, mlist.begin() + 55, 3);     // n in [18,55], m = 3
    std::fill(mlist.begin() + 55, mlist.begin() + 59, 4);     // n in [56,59], m = 4
    std::fill(mlist.begin() + 59, mlist.begin() + 76, 3);     // n in [60,76], m = 3
    std::fill(mlist.begin() + 76, mlist.begin() + 239, 4);    // n in [77,239], m = 4
    std::fill(mlist.begin() + 239, mlist.begin() + 247, 5);   // n in [240,247], m = 5
    std::fill(mlist.begin() + 247, mlist.begin() + 284, 4);   // n in [248,284], m = 4
    std::fill(mlist.begin() + 284, mlist.begin() + 991, 5);   // n in [285,991], m = 5
    std::fill(mlist.begin() + 991, mlist.begin() + 1007, 6);  // n in [992,1007], m = 6
    std::fill(mlist.begin() + 1007, mlist.begin() + 1083, 5); // n in [1008,1083], m = 5
    std::fill(mlist.begin() + 1083, mlist.begin() + 2015, 6); // n in [1084,2015], m = 6
    std::fill(mlist.begin() + 2015, mlist.begin() + 2031, 7); // n in [2016,2031], m = 7
    std::fill(mlist.begin() + 2031, mlist.end(), 6);          // n in [2032,2204], m = 6

    return mlist;
}

/*
 * @brief 使用长除法算法进行Chebyshev多项式除法
 *
 * @param f 被除数（多项式f的系数）
 * @param g 除数（多项式g的系数）
 * @return longDiv 返回长除法结果，包括商和余数
 */
longDiv Bootstrapper::LongDivisionChebyshev(const vector<double> &f, const vector<double> &g)
{
    uint32_t n = Degree(f);
    uint32_t k1 = Degree(g);

    if (n != f.size() - 1) {
        cout << ("LongDivisionChebyshev: The dominant coefficient of the divident is zero.");
    }

    if (k1 != g.size() - 1) {
        cout << ("LongDivisionChebyshev: The dominant coefficient of the divisor is zero.");
    }

    if (int32_t(n - k1) < 0)
        return longDiv(std::vector<double>(1), f);

    std::vector<double> q(n - k1 + 1);
    std::vector<double> r(f);
    std::vector<double> d;
    d.reserve(g.size() + n);

    while (int32_t(n - k1) > 0) {
        d.clear();
        d.resize(n + 1);

        q[n - k1] = 2 * r.back();
        if (IsNotEqualOne(g[k1]))
            q[n - k1] /= g.back();

        if (int32_t(k1) == int32_t(n - k1)) {
            d.front() = 2 * g[n - k1];

            for (uint32_t i = 1; i < 2 * k1 + 1; i++) {
                d[i] = g[std::abs(int32_t(n - k1 - i))];
            }
        } else {
            if (int32_t(k1) > int32_t(n - k1)) {
                d.front() = 2 * g[n - k1];
                for (uint32_t i = 1; i < k1 - (n - k1) + 1; i++) {
                    d[i] = g[std::abs(int32_t(n - k1 - i))] + g[int32_t(n - k1 + i)];
                }
                for (uint32_t i = k1 - (n - k1) + 1; i < n + 1; i++) {
                    d[i] = g[std::abs(int32_t(i - n + k1))];
                }
            } else {
                d[n - k1] = g.front();
                for (uint32_t i = n - 2 * k1; i < n + 1; i++) {
                    if (i != n - k1) {
                        d[i] = g[std::abs(int32_t(i - n + k1))];
                    }
                }
            }
        }

        if (IsNotEqualOne(r.back())) {
            // d *= f[n]
            std::transform(d.begin(), d.end(), d.begin(),
                           std::bind(std::multiplies<double>(), std::placeholders::_1, r.back()));
        }
        if (IsNotEqualOne(g.back())) {
            // d /= g[k1]
            std::transform(d.begin(), d.end(), d.begin(),
                           std::bind(std::divides<double>(), std::placeholders::_1, g.back()));
        }
        // f-=d
        std::transform(r.begin(), r.end(), d.begin(), r.begin(), std::minus<double>());
        if (r.size() > 1) {
            n = Degree(r);
            r.resize(n + 1);
        }
    }

    if (n == k1) {
        d = g;

        q.front() = r.back();
        if (IsNotEqualOne(g.back())) {
            q.front() /= g.back(); // q[0] /= g[k1]
        }
        if (IsNotEqualOne(r.back())) {
            // d *= f[n]
            std::transform(d.begin(), d.end(), d.begin(),
                           std::bind(std::multiplies<double>(), std::placeholders::_1, r.back()));
        }
        if (IsNotEqualOne(g.back())) {
            // d /= g[k1]
            std::transform(d.begin(), d.end(), d.begin(),
                           std::bind(std::divides<double>(), std::placeholders::_1, g.back()));
        }
        // f-=d
        std::transform(r.begin(), r.end(), d.begin(), r.begin(), std::minus<double>());
        if (r.size() > 1) {
            n = Degree(r);
            r.resize(n + 1);
        }
    }
    q.front() *= 2; // Because we want to have [c0] in the last spot, not [c0/2]

    return longDiv(q, r);
}

/*
 * @brief 计算多个密文加权和
 *
 * @param ciphertexts 输入的密文向量
 * @param constants 加权系数
 * @return Ciphertext 返回加权和的密文
 */
Ciphertext Bootstrapper::EvalLinearWSumMutable(vector<Ciphertext> &ciphertexts, const vector<double> &constants)
{
    long minLevel = ciphertexts[0].l;
    int minIndex = 0;
    int wssb = 1;
    if (wssb != 0) {
        for (int i = 1; i < ciphertexts.size(); i++) {
            if (ciphertexts[i].l < minLevel) {
                minLevel = ciphertexts[i].l;
                minIndex = i;
            }
        }
    }
    for (int i = 0; i < ciphertexts.size(); i++) {
        AdjustLevelsAndDepthInPlace(ciphertexts[i], ciphertexts[minIndex]);
    }

    Ciphertext weightedSum = scheme.MultByConst(ciphertexts[0], constants[0]);
    Ciphertext tmp;
    for (uint32_t i = 1; i < ciphertexts.size(); i++) {
        tmp = scheme.MultByConst(ciphertexts[i], constants[i]);
        scheme.AddAndEqual(weightedSum, tmp);
    }
    scheme.ReScaleByAndEqual(weightedSum, 1);

    return weightedSum;
}

/*
 * @brief 调整两个密文的级别和深度，使其一致
 *
 * @param ciphertext1 第一个密文
 * @param ciphertext2 第二个密文
 */
void Bootstrapper::AdjustLevelsAndDepthInPlace(Ciphertext &ciphertext1, Ciphertext &ciphertext2)
{
    if (ciphertext1.l < ciphertext2.l) {
        scheme.ModDownByAndEqual(ciphertext2, ciphertext2.l - ciphertext1.l);
        ciphertext2.scale = ciphertext1.scale;
    } else if (ciphertext1.l > ciphertext2.l) {
        scheme.ModDownByAndEqual(ciphertext1, ciphertext1.l - ciphertext2.l);
        ciphertext1.scale = ciphertext2.scale;
    } else {
        return;
    }
}

/*
 * @brief 应用双角迭代算法对密文进行处理
 *
 * @param ciphertext 输入的密文
 * @param numIter 迭代次数
 */
void Bootstrapper::ApplyDoubleAngleIterations(Ciphertext &ciphertext, uint32_t numIter)
{
    int32_t r = numIter;
    for (int32_t j = 1; j < r + 1; j++) {
        scheme.SquareAndEqual(ciphertext);
        scheme.ReScaleByAndEqual(ciphertext, 1);
        scheme.AddAndEqual(ciphertext, ciphertext);
        double scalar = -1.0 / std::pow((2.0 * M_PI), std::pow(2.0, j - r));
        scheme.AddConstAndEqual(ciphertext, scalar);
    }
}


void Bootstrapper::BootStrapping(Ciphertext& cipher)
{
    // 10.对密文进行缩放
    double smaller = 0.00195313;
    if (context.N == 32768) {
        smaller = 0.015625;
    }
    Ciphertext ct = scheme.MultByConst(cipher, smaller); // 最开始除以512进行缩放
    scheme.ReScaleByAndEqual(ct, 1);

    // 11.模升
    scheme.ModUpToQAndEqual(ct);

    // 12. CoeffsToSlots

    double pre = 0.5;
    double k = 1;
    double N = 1 << context.logN;
    double constantEvalMult = pre * (1.0 / (k * N));

    ct = scheme.MultByConst(ct, constantEvalMult);
    scheme.ReScaleByAndEqual(ct, 1);              

    Ciphertext ctxtEnc = encodingMatrix.EvalCoeffsToSlots(encodingMatrix.mU0hatTPreFFT, ct); 

    Ciphertext conj = Ciphertext(ctxtEnc.N, ctxtEnc.slots, ctxtEnc.l, ctxtEnc.scale);
    scheme.Conjugate(conj, ctxtEnc);

    Ciphertext ctxtEncI = scheme.Sub(ctxtEnc, conj);

    ctxtEncI = scheme.DivByi(ctxtEncI);

    ctxtEnc = scheme.Add(ctxtEnc, conj);

    // step 4------------EvalMod------------
    std::vector<double> coefficients; 
    coefficients = gCoefficientsSparse;
    Ciphertext ctxtEnc2 = EvalChebyshevSeriesPS(ctxtEnc, coefficients, -1, 1);

    ApplyDoubleAngleIterations(ctxtEnc2, 3);


    Ciphertext ctxtEncI2 = EvalChebyshevSeriesPS(ctxtEncI, coefficients, -1, 1);
    ApplyDoubleAngleIterations(ctxtEncI2, 3);

    // step 5----------SlotsToCoeffs--------

    ctxtEncI2 = scheme.MultByi(ctxtEncI2);

    Ciphertext CtSTC = scheme.Add(ctxtEnc2, ctxtEncI2);

    scheme.MulScalaInt(CtSTC, 2);

    CtSTC = encodingMatrix.EvalSlotsToCoeffs(encodingMatrix.mU0PreFFT, CtSTC); // 密文L = 15-->12

    scheme.ReScaleByAndEqual(CtSTC, 1);

    double bigger = 512;
    if (context.N == 32768)
    {
        bigger = 64;
    }

    scheme.MulScalaInt(CtSTC, bigger); // 对密文进行缩放
    cipher = CtSTC;
}